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SOME BASI C MATHEMATI CAL METHODS 
OF DIFFUSION THEORY 

I. INTRODUCTION 


The purpose of this report is to provide an Introductory treatment of the 
basic mathematical aspects of diffusion theory as It applies to the atmosphere, 
starting with molecular diffusion and leading up to the statistical methods of 
turbulent diffusion. The scope of the treatment is evident in that hydrodynamlcal 
theory, which is a necessary foundation of diffusion theory when applied to a 
medium as complex as the atmosphere, has been kept to a minimum. The basic 
concepts and equations of diffusion are developed on an elementary level, with 
emphasis on the mathematical steps. Many of the derivations given cannot 
easily, if at all, be found in the literature on diffusion. The report is not meant 
to be definitive in any sense. It is hoped that sufficient material has been 
treated to provide a background for further reading to understand today' s 
specialized papers in the field. 
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II. BASIC EQUATIONS OF DIFFUSION 




A. Stationary Medium 

Diffusion is the process whereby a substance introduced in a localized 
region of some medium (e.g. « a fluid) spreads throughout the medium either 
by molecular motion or by turbulence. The basic equations of molecular dif- 
fusion are derived first. Turbulent flow will be treated later. If x denotes the 
concentration of a diffusing substance, defined as the quantity of substance per 
unit volume, then relative to a set of rectangular axes the vector rate of 
transfer of substance through unit area of a section of medium is called the 
flux F . In many diffusion problems the assumption made is that F is pro- 
portional to the concentration gradient normal to the section. That is. 


F = 

X 



F =-D 

y 


ay 


F = 
z 


-D 


dz 


or in vector form 


F= -D^X 


( 1 ) 


where D = D(x,y,z,t,x) Is the diffusivity or coefficient of diffusion which may, 
in general, be a function of the coordinates, the time, and the concentration. 
The quantity V x is the gradient of the concentration defined by 


Vx = 


3x 


ay 



k 


f 


where i , j , and k are unit vectors along the x, y, and z axes, respectively. 

The negative sign indicates that the diffusion occurs in the direction of decreasing 
concentration. 
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r 
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We derive next the equation connecting the time rate of change of the 
concentration with the spatial rate of change of the concentration at any given 
point. This is accomplished by considering an arbitrary volume of fluid V, 
bounded by the surface S, as shown in Figure 1, and applying the principle of 
continuity to the diffusing substance. 


z 



Figure 1. The unit normal and flux vectors at a point 
on a surface S bounding the volume V. 


Let M be the total mass of diffusing substance within V. It is clear that 
M will be a function of time. From the definition of x it follows that 


M(t) = /^x dV 


( 2 ) 


3 
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At every point of the surface S there is a flux F making some angle 
with the unit normal *n at the point. The only flux out of the volume V has to 
be along the direction of n, which is just the component of F parallel to n, or 
F • n. Therefore, the total flux leaving the volume V is obtained by integration 
over the surface S. The total flux has units of mass per unit time and, there- 
fore, may be called the current, which is denoted by I(t) . Thus 


This surface integriil may be converted to a volume integral by use of the 
divergence, or Gauss’ , theorem which states that 


f_ F • *n ds = Lv • F dV 
■ D ' V 


and so 


I(t) = f^V • F dV 


Since I is a mass current, we have 




which, when combined with equations (2) and (3), results in 


/„V.?dV=-|f /yXClV 


V 


or 




dV= 0 


Since the volume is arbitrary, the integrand must vanish, yielding 


8l 


-V . F 


(4) 


which Is simply the expression for the conservation of diffusing matter, namely 
that the rate of change of the concentration at any poii.t must equal the net flux 
at the point. The assumption of conservation is valid provided there is no 
creation or destruction of the diffusing matter. 

In rectangular coordinates, equation (4) assumes the form 


„ 3F 8F 8F 

2L . — 2. 

3t dx 9y 9z 


(5) 


Substituting equation (l) into equation (4) yields the relation we sought con- 
necting the time rate of change of x with its spatial rate of change. 


9t 


= V 


D^X 


( 6 ) 


which in rectangular coordinates is 


_ d 
9t 9x 




(7) 
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Equation (6) is the fundamental equation governing the diffusion process in 
stationary isotropic media. If D is constant, it reduces to 


. ( 8 ) 


where V^x ^ V • Vx is the Laplacian of x which, in rectangular coordinates, has 
the form 



( 9 ) 


When D is constant, the diffusion, as expressed by equation ( 8) , is called 
Fickian diffusion. The formulation given here is strictly phenomenological, 
resting on the assumption that the flux is proportional to the concentration 
gradient as expressed by equation (l) . 


B. Moving Medium 

In the case of a fluid moving with a laminar flow velocity u, as might be 
the case in the atmosphere when there is a wind, the flow will contribute to the 
flux by physically transporting the diffusing matter. The additional flux is 
iTx, and the flux equation (l) becomes 


F 


= ux - X 


( 10 ) 


Substitution of equation ( 10) into equation ( 4) gives 


If = -V • (ITx -D^X) 

= -V • ( ux) + V • (D^x) 


If u ts constant, the second term is zero, which implies incompressible flow. 
Taking D to be constant also, we have 


at 


+ u • Vx = DV*x 


( 11 ) 


which is the fundamental equation of diffusion commonly used in the case of a 
fluid w th laminar flow. 

In a situation where D is a function of time only, D « D(t) , the problem 
can be handled by defining a new time scale t by 


dT = D dt 


Then, by differentiation. 


^ a ^ ^ ^ £) 

at 3t dt 8t 


_ ^ dx dr 
\ dt dr dt \ 


Similarly, 


u=vD , u=vD 
y y z z 


so that equation (ll) becomes 





which Is a form identical to equation (ll) except that D does not enter the 
equation. It should be clear that *v is the velocity in terms of the new time 
scale T. For a known u and D( t) « the value of *v is given by 



C. Anisotropic Medium 

In all of the equations developed thus far» the diffusion coefficient has 
been assumed to be the same in all directions. When this is not the case» the 
medium is said to be anisotropic. The equations for the flux as given by equa- 
tion (lO) in the general case of a moving medium must now be modified. It is 
convenient at this point to express the flux equations (10) and succeeding equa- 
tions in rectangular coordinates. By a coordinate transformation the equations 
may then be cast in terms of any other system of coordinates. 

The new equations for the flux are 


F 

X 


F 

y 


F 

z 


u D ^-D ^-D ^ 

X* 11 ax 12 ay 13 ax 


„ D ^-D ^-D ^ 

y^ 21 dx 22 8y 23 dz 


ux-D -^-D -^-D ^ 

tT 31 3x 32 8y 33 9z 


f 


which can be written concisely in tensor notation as 


F = u V - 

i r 


) ^ 
ij 


(13) 
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with summation implied on the repeated subscript and with the understanding 

that F , F , F are replaced with F • F , F , respectively, and x, y, z are 
X y z 1 4b u 

replaced with x , x , x . 

1 / u 

The quantity Is the tensor diffusivlty and consists of nine quantities, 
which can be written in the matrix form 



D 

D, 

D 


1 

11 

D 

D „ 

12 

13 

1 

21 

D 

D 

22 

23 

> 

D 

D 

31 

32 

33 


(14) 


It is evident from equation (13) that the diffusion in a particular direction 
depends not only on d\/ 9x but also on 8x/ 9y and 9x/ 9z. 

Substituting equation (13) into equation (4), 


9t 


-V 


9F. 


9x 


i 


_9 

9x 




9u 9D V 

9Xj I 9x, 8Xj 8x^ ij 9Xj8x^ 


(15) 


9 






Assuming and are constant, equation ( 15) reduces to 


_9!i_ 

at i ax. ij ax. ax. 

i 1 j 


(16) 


Because of the difficulty in finding solutions of the general equation (15), v:: 
shall restrict ourselves to the more useful particular form, equation (16) . By 
transformation of coordinates to a new set of rectangular axes, 1 1 , ^ 2 » ^3 
( called principal axes) , equation ( 16) can be simplified to 


9l_ 

at 


= -U. 


- A 

ait. i 




(17) 


or, in expanded form. 


„ iiL ,, lx. 
at ‘ail" ^9^2 



. (18) 


Evidently, the transformation has eliminated the cross derivatives and reduced 
the diffusion iensor to a diagonal form: 



0 


D.. = 
ij 


0 


0 



(19) 


The values of the coefficients Dj, D 2 , D 3 (known as principal diffusion coeffi- 
cients) in general will depend on the former coefficients D„ given in equation 

(14) and also on the coordinate transformation coefficients. 
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We might note at this point that the equations developed herein are 
linear differential equations if D ( or the in the case of anisotropic media) 

are constants or functions only of the coordinates and time. The property of 
linearity is important because it Implies the principle of superposition; i.e. , 
if one or more solutions are found the sum of these solutions is also a solution. 
It is not at all evident that the boundary conditions associated with a given 
problem can be satisfied. In some cases a particular solution will be sufficient, 
but generally a sum, usually an infinite sum, of particular solutions will be 
necessary to satisfy bovindary conditions. 

In summary, we list here the equations most tractable to solution in a 
practical situation together with the conditions under which they are applicable: 


Equation ( 8 ) : 


9t ^ 


Isotropic media, constant O 


Equation (ll); 

Isotropic media, constant D, laminar 
flow with constant u 


Anisotropic media, constant D., 
laminar flow with constant u 

In this last equation we have reverted to the standard symbols Xj, X 2 , X 3 to 
denote the coordinates, with the understanding that these are principal 
coordinates. 


+ u • Vx = DV^x 


Equation (17) : 


dx 
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III. SOURCE SOLUTIONS 


A. Instantaneous Sources 

The diffusion equations ( 8) , ( 11) , and ( 17) can be solved by the standard 
method of separation of variables or by the Laplace transformation. The litera- 
ture on this is quite extensive [1-3] . The purpose here is to consider some 
solutions which are relevant to the atmosphere and which form a guide for the 
statistical theory of turbulent diffusion. The simplest example is the one- 
dimensional problem of a medium at rest (u = O) . If the diffusion is along one 
direction only, this implies that the concentration is imiform in the other two 
directions, meming that the concentration gradient everywhere along these two 
directions is zero. 

Denoting the diffusion direction by x, one can easily verify by direct 
substitution that the following expression, which can be derived by use of the 
Laplace transformation, is a solution of equation ( 8) ; 


X = 


-sTt 



( 20 ) 


where A is an arbitrary constant. If we integrate this from -® to +® , we 
obtain the total amount Q of diffusing substance per unit cross-sectional area 
of a cylinder of infinite length and cross section: 


CO CO - ■ 

Q= f Xdx - f e ^^^dx = 2 An/ ttD 


From this it follows that 
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T 
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therefore, equation (20) becomes 




k 


; / 
/ 


I 



( 21 ) 


This result describes the diffusion of a quantity of substance Q which was 
initially concentrated in the plane x = 0, known as an instantaneous plane source. 
The fact that this is an idealized source is evident by observing that y » as 
t-* 0. Nevertheless, this solution, as we shall see, forms the basis for obtain- 
ing more realistic solutions. 

Since the solution, equation (2l), has the functional form of a Gaussian 
or normal curve, we may compute the second moment of x, defined by the 
integral 


^ ^ dx= 2QDt 

-OD -CO 


This result characterizes the "spread" of the concentration x about its maximum 
value, which occurs at the origin. If we take this result and divide it by the 
total amoimt of diffusing matter Q, we obtain a quantity a defined by 


o-*= 2Dt 


( 22 ) 


The quantity a has the dimensions of a length and defines a kind of root-mean- 
square distance to which the substance has diffused. To make this clear, a 
plot of equation (21) for several values of <r^ is shown in Figure 2. It is evident 
that the spread increases as a increases. This behavior justifies the use of o 
as a parameter characterizing the spread of x. In statistical work, a serves 
as a convenient scale of the width of the distribution and is known as "standard 
deviation. " Its square is called the variance. Note that or is a function of the 




t 



C 
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dlffusivity and time. In terms of a, equation (21) is 


X 


Q 



( 23 ) 


As pointed out earlier, the instantaneous plane source is an idealized 
construct. A more realistic approach is to consider a source distributed over 
a finite region of x. 

Such a formulation would be 


X = Xo for - a < X < a , and t = 0 
X = 0 elsewhere, t = 0 


where 2a is the thickness of the cloud. Uniformity along the y and z axes is 
again assumed. If we divide the interval - a < x < a into an infinite number of 
infinitesimally thin sheets, each sheet can be considered an instantaneous plane 
source as defined previously. Then, because the diffusion equation is linear, 
the solution is obtained by summing up all the plane-source solutions arising 
from each sheet. 

If x* denotes the coordinate of one of these thin sheets, the distance of 
a point X from the sheet is x - x* , where dx' is the differential thickness of the 
sheet. Thus, each sheet contributes an amount dx to the total x< That is, 

(x-xM^ 
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Integrating from -a to +a, we obtain the solution 


I 



(x-x’)2 

r 2ff^ 

.1 ® 


( 24 ) 


It is possible to write this integral in terms of the error function if we define a 
new variable by 




t 


and, denoting the new limits on the Integral by p and q, 



( 25 ) 
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f 
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where erf p and erf q represent the error functions of p and q as defined by the 
preceding integrals. Tables exist for evaluating the error function for a given 
value of its argument. Moreover, the error function possesses the following 
properties: 


erf ( “p) = -erf p , erf (O) = 0 , erf(«>) = l 


Figure 3 presents a plot of equation ( 26) for several values of the parameter 
Dt/ a}. At this point it might be well to remember that the time dependence of 
X is contained in the quantity cr defined by equation (22) . 

The maximum value of x at any given time occurs at x = 0, the center of 
the cloud. With x = 0 we have 


a 



17 


t 



I 


i 



x-a 
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therefore, equation (26) yields 


Y = Y„ erf 
^max 0 


(^/p) 


As an example of the use of the foregoing formulas, let the imlf-thlckness 
a of the cloud, the dlffuslvity D, and the initial concentration be given by 


2 1 - 

a=lm , D=0.66— , y =0.60-^ 

sec 0 m 


and suppose we seek the concentration at a distance x = 60 m when t = 1800 sec 
(30 min) . We first compute o. 


(T = */2Dt = '7210. 66 ~ 1(1800 sec) = 48.7 m , 

\ S0C / 


then p and q. 


p, ^2-!4- ,0.89 ., = -12^=0.86 

■72 (48.7) «J2(48.7; 


and finally, with the aid of a table of the error function, equation (26) gives 


= (0.792 - 0.776) = 0.0048 


2 m 


V 




\ 





k 


In addition, one might compute also the maximum concentration of the distribu- 
tion for this given time of 1800 sec. From equation (27) this is 


= = 0.010 

max m m 


Note; Any consistent set of units may be used. For example, if 
length is measured in meters and time in seconds, D must be 
in m^/sec. K length is measured in centimeters, then D should 
be in cm*/ sec. Similarly, the concentration may be in kg/m* 
or gm/ cm*. A larger unit for the time may be used provided D 
is expressed in this larger unit. For example, if D = 0. 2 
cm*/ sec, it would convert to D = 0. 2 cm*/ sec x 60 sec/ min, or 
D = 12 cm*/ min. 

Thus far we have treated the problem of a source extended in a finite 
region in an infinite mediiim. For this problem, therefore, the solution tends 
to zero as x approaches infinity. If, however, the medium is not infinite, then 
we must assume the existence of a boundary at some distance x = L. At such 
a boundary, we can have either total reflection or total absorption of the dif- 
fusing matter or, more generally, some reflection and some absorption, all 
depending on the nature of the boundary. In the case of the atmosphere, this 
boundary may be the ground, a hillside, a building in the path of the diffusing 
matter, or it may even be another layer of the atmosphere which possesses 
different diffusion characteristics. 

If at the boundary we have a prescribed flux, the mathematical statement 
of the flux should be a linear function of x* If not, a solution can rarely, if 
ever, be foimd. An example of a linear flux condition is 


-D 


an 


= H (x - Xo) at x = L 


f 


where the left member is the definition of flux (n specifies the normal direction 
at the boundary surface) . The right side, where II is some constant, states that 
the flux is proportional to the difference in concentration between the surface 



T'— “ 

I 


1 
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and the surrounding medium at constant concentration xo« case of heat 

flow* this condition is an expression of Newton* s law of cooling* valid for small 
temperature differences. 

If* however* the flux is totally reflected at the boundary surface* we 
have the condition of zero flux: 

D-?^ = 0atx=L 
dn 


A solution with this boundary condition can be easily determined by the 
method of images. For example* we imagine that there is a source at x = 2L 
which is a mirror image of the real source, i.e. identical in every respect. 
The two sources diffuse in identical fashion* and midway between at x = L at 
the boundary surface* the flux from one exactly cancels the flux from the other* 
resulting in the prescribed zero flux. The geometry of this is depicted in 
Figure 4. Furthermore* since the diffusion eqxiation is linear* the sum of 
the two solutions will yield the desired solution to our boundary-value problem 
with the impermeable boun-'ury. 



-a 0 +a 
SOURCE 


L 


X 


2L-a 2L 2L+a 
IMAGE 


Figure 4. The real source and its image, where x = L is an impermeable 
boundary and ail quantities are measured relative to a coordinate x 
whose origin is at the center of the source. 
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The solution of the mlrror-Itnage source will be an Integral similar to 
equation (24) for the actual source. That integral was obtained by placing the 
origin of the coordinate system at the center of the source. We do the same 
with the image source except now the new coordinate is denoted by x. Thus, 
this solution is 


X = 


(X-X’)2 

^ x=a - - V- 

XO j- ^ 2(7“ 

2ira^ jc=-a 


dx» 


The coordinate transformation connecting x and the original coordinate x is 


X = X - 2L 


With this relation, the integral can be written relative to the original coordinate. 
The relevant quantities transform as follows: 


X - x’ = (x - 2L) - (x* - 2L) 
=* X - x’ 
dx’ = dx , 


and the lower and upper limits become 2L - a and 2L + a. Thus, the mirror- 
image solution is 


Xo 

•J 2ita^ 


2 1/+*a 
/ 


2L-a 
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This Integral is identical to equation (24) except for the limits. Proceeding 
as before, we find 


pt = 


X + a - 2L 


= 


X - a - 2L 

4 2(t^ 


and the solution Is 


X = “ {erf p’ - erf q’ } 


Finally, the complete solution la the sum of this and equation (24) ; 


Xo 

X = ” { erf p - erf q + erf p' - erf q' } 

Ct 


or, In explicit form, 


X = 


Xo 


,x + a ,x-a ,x + a-2L ,x-a-2L 
erf erf + erf — - erf 


Vi? n/J? 


^/2?’ n/2? 


. (28) 


At the boundary x = L, 


, - V , . a + L , a - L 

X (x = L) = Xo ) erf IT + erf r: 

42a^ 42a^ 


(29) 
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We consider next two- and three-dimensional problems. It can be 
verified that a source solution of the two-dimensional diffusion equation. 



w'here again we assume the flow velocity to be zero, is given by 


X = Ae 



where A is an arbitrary constant. For the diffusion to be indeed two-dimensional 
everywhere, it must be presumed that there is uniformity of concentration in the 
z direction and that the source is uniformly distributed over the entire z axis. 
Therefore, the constant A can be determined by demanding that the total amount 
of diffusing matter, Q^, in a cylinder of finite height concentric with the z axis, 
but of infinite radius, must be conserved. This is given by 


/ 


Z2 ® ® 

Qj ^ / / X dx dy dz 

2 j -00 -00 


X 


^2 00 ^ 2 ® ** 

= A J dz f e ^ dx f e 

Zi -CO -CO 


o»2 


dy 


= A (z2 - Zj) ( 42na^ ^2ir<x^) 
- 2ir(7^A (z2 - zi) 
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i ^ = 27r<T*A 

I Z2 - zi 

ii 
J 

i 

; . Since zj - Z| is the height of the cylinder, the quantity on the right is the total 

^ matter per unit length of cylinder which, to be conserved, must equal the initial 

line density Q on the z axis. Thus 


Q = 2ir(T^A 


or 


A = ^ 

2it(T^ 


so that 



(30) 


Note the dimensions of Q: the amount of matter per unit length 
of line. 

The solution, equation ( 30) , governs the diffusion of an instantaneous and 
infinitesimally thin source located on the z axis; hence, the name **line source” 
is given to it. Converting to the polar coordinate p defined by 


i 

i 

t 

I 




\ 
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the solution can be written as 


I 


X 




(31) 


which, with the exception of the constant Q/ 27r(r^ (constant with respect to the 
coordinates, since a is still a function of time) , has exactly the same form as 
the one-dimensional solution. The role of x is now played by p. Thus, at any 
given time the concentration at points equidistant from the origin is the same; 
i. e. , the contours of constant x are circles centered at the origin. Any plane 
parallel to the z axis and through the origin would show a concentration profile 
exactly as in the one-dimensional case. 

The solution for a long cylindrical source of radius a cannot be obtained 
in terms of elementarj' functions. We merely quote the result which is 


_JSL 

Xo 2(t^ r 2o^ _ / \ j , /no\ 

X = p-e I e ^(“^jp’dp’ (32) 


where p' represents the radial coordinate of a point in the cylinder (the inte- 
gration variable over the region) , p is the observation point, and Iq is the modified 
Bessel function of the first kind of zero order. There is no way of evaluating 
this integral except numerically. An elementary solution, however, exists for 
points on the axis of the cylinder. At such points p = 0 and 1^ (O) = 1. By 
elementary integration, we obtain 


X = 




(33) 


This simple formula is useful in determining the decay of the concentration on 
the axis, where the concentration is always above that of surrounding points. 
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In three dimensions the instantaneous point source solution it: 


X = Ae 




In this case the integral over all space must equal the initial source strength, 
which we now denote by xo instead of Q because it has the same dimensions as 
X. Thus 


00 00 GO 

X. = A/ f f 

mGO —00 —00 





dx dy dz 


.A / 

-00 



z 


2 



e dz 


= A(2jto-*)^^ 


or 


A = 


Xo 


(2ircr^) 


i,V! 


therefore. 


x = 


Xo 


+ z^ 
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( 34 ) 
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where Is the radial coordinate. Again we see that the cloud 

grows along any diameter, as in the one-dimensional case. 

The point source solution, equation ( 34) , may be Integrated to obtain 
the solution for a uniform spherical source of radius a. In spherical coordinates 
with axes at the center of the sphere, we let r be the radial coordinate to the 
point of observation, r* the coordinate of a volume element of the sphere, and 
R the distance to the observation point. Because the concentration Is spheri- 
cally S3'mmetric, without loss of generality we cu^se the observation point on 
the z axis. 

Referring to Figure 5, we see that 


r" = r^ + r ♦ ^ - 2rr* cos 0 


1 



Figure 5. The geometry for the integration of the point source 
solution for a spherical source of radius a. 
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We now make a change of variable to u in the first integral and to v in the 
second: 


r - r’ r + r’ 

^ v = 

•/2 cr ^ a 


Writing also 


r + a r - a 

p = q = 

*7T ff cr 
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we have 


Xo 


X = 


-J 


(r -»/2’ o-u) e du 


i r 

n/T(T 


P -V* 

+ f (r -*'/T a v) e dv 

r 

^(T 


These integrands are identical* u and v being variables of integration. We can* 
therefore* interchange the limits on the first integral and then write the sum of 
the two integrals as a single integral* which is allowed by the fundamental 
property of definite integrals for continuous intervals. Thus 


j ( r - n/T (t u) e du 

q 


r -u I — r ~u 

r I e du - n/ 2 (T I ue du 

q q 

We now write the first integral as the difference of two integrals* one from 0 
to p and the other from 0 to q. The second integral is integrated directly by 
elementary means. Thus 



Xo 
Vir r 


Xo 

n/1t r 
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(35) 


which is the solution for a sphere of radius a. 


B. Continuous Sources 

We consider now the three fundamental sources, plane, line, and point, 
when they emit continuously for times t > 0 at the rate of Q(t’ ) units per unit 
time. In the time element from t’ to t’ + dt’ , an amount Q(t’ ) dt’ of material 
is emitted. Moreover, the standard deviation of the distribution will depend on 
the time t’ in the following way; 


(T^ = 2D(t -t») 


(36) 


From equation (23) the contribution to x due to an elemental plane 
source Q dt’ is 


Q dt’ e 

dx = 


4D(t-t’) 


2 '^ ( t - t ’)‘/2 


Integrating this over the time of emission, we have 



4D(t-t’) 

(t-t ’)‘/2 


dt’ 
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which is the general solution if the emission rate Q is some general function of 
time. For a constant emission rate, this becomes 


00 -s 
e 


= ^ f ^ds 
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where we put 


s = 


2n/ D(t - t» ) 
Integrating by parts. 
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where erfc a is the complementary error function defined by 


32 


V 


T" 


T 

I 


« ® 9 

- 2 r -8^ 

erfc Of = J e ds 
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= 1 - erf a 


since the first integral evaluates to unity and the second is the error function, 
by definition. 


Similarly, employing equation (31), the solution for a continuous line 
source is 


X = — fQ(V)~ 

^ 47tD •' ' t - 


4D(t-t» ) 


which for constant Q is 


® -s 

-7^ f — is 

4ttD •{ s 


where in this case we defined s by 


4D(t - t’) 
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The exponential integral involved in this soluticm cannot be evaluated in closed 
form. Its series expansion is 
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where y - 0 57721 .... is Euler* s constant. Thus» 
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For large values of t the series term can be neglected, giving the approximate 
form 


X ' 




47tD 


(.„ ^ - y) 


(39) 


Finally, in like fashion, we can write the solution for a continuous point 
source with the aid of equation ( 34) : 


X = 


t ^ 4D(t-t») 

f Q(t*) r~ dt» 


8 ( 710 )^-^^ 0 


(t-f) 


V2 


I’uttlng 
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and for constant Q, we have 


Dr r 
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= - 2 - 
47rDr 
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We might note for this solution that as t -* « it reduces to Q/ 4irDr which is a 
steady state distribution. 


C. Source Solutions in a Moving Medium 

The extension of the source solutions to the case of a moving medium is 
rather simple. We note that the solutions obtained in the case of a stationary 
medium may be considered solutions as observed from a coordinate system 
moving with the medium. Then, by transforming coordinates to a stationary 
system, the corresponding solutions for a medium in motion relative to the 
system at rest are obtained immediately. 

If u, V, w are the vector components, assumed constant, of the velocity 
of the medium, and x’ , y’ , z’ are the coordinates of a point as observed in the 
stationary system, then the coordinates of the point relative to a moving system 
are given by 


X = X* - ut 


y = y» - vt 


z = z’ - wt 





Substitution of these expressions into the equations for instantaneous sources 
for a stationary medium will yield the solutions for a moving medium. For 
example, the point-source solution, equation (34), would take the form 


XO 

(2jr)’'^*or® 


(x*-u t)^_ (y»-vt)^. (z»-wt)^ 




2<t* 
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with similar results for the plane and line sources. Ckie can always choose a 
coordinate system such that one of the axes, e.g. the X axis, is parallel to the 
velocity. For the point source we would have the simpler result: 


X = 


Xo 


2a^ 2a* 2a* 


(42) 


where the prime on x may be dropped with the understanding that it is measured 
from a fixed system. The fact that these expressions satisfy the diffusion equa- 
tion for a moving medium can be verified by direct substitution. 

In general, if x (x* , y* , z’ , t) is any instantaneous source solution of the 
diffusion equation for a stationary medium 




where the prime on the Laplacian indicates derivatives with respect to the primed 
coordinates, then y(x, y, z, t) is a solution of the diffusion equation for a moving 
medium 


ll + V . 

at 


VX = 


where x, y, z are given by the relations of equation (41) and V = (u, v, w) . 



From the solution of equation (42), the expression tor a continuous point 
source in a wind can be determined. This would represent roughly a model of a 
chimney plume, because the plume may be considered as a continuous series of 
point source clouds or puffs. At time t* the source emirS an elemental amount 
Qdt* of matter, while an element of air at the point ( x, y, z) at time t will have 
been initially at x ~ u(t - t*) » y* z because of the wind. Thus, the concentration 
dji at (x, y, z) and at time t due to an instantaneous puff of amount Qdt' is 


dx = 


Qdt* 
8(jrD(t-t*)] 


[x-u(t-t* ) 1^+y^+z^ 
4D(t-t») 


V2 


If the emission rate Q is constant, the total concentration will be 


X 


t . tx-u(t-t*)l V4-z^ 

r 1 df 


For a source maintained indefinitely, which is most important in practice, the 
limits on this integral would be from 0 to » . hi this case the integral can be 
evaluated in closed lorm. Thus, 


X = 


ux 





00 


/ 



uV 

16D^s* 


) 


ds 



47rDr ® 


(43) 
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where r* = x* + y* + z* and 


s 


r 

2*/D(t-t») 


We might note that equation (43) is independent of time, as expected for 
a source maintained indefinitely at a constant emission rate. 

Furthermore, studies on smoke clouds have shown that the cloud has the 
form of a long, thin plume if the wind velocity is not too low. In this case one 
is usually interested in the concentration values near the axis of the plume, 
where y = z - 0. On the axis itself the above expression takes the very simple 
form 


^ 47TDX 


For points near the axis relatively far downwind the quantity (y^ + z^)/x* is 
small. With the aid of the binomial expansion and neglecting powers of this 
quantity higher than the first, 



Putting this into the exponent of equation (43) gives the approximate form 


Q 

X = TfiC 


u(v^->-z^) 

4Dx 
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One can also, in practice, replace the r in the coefHcient of the exponential 
factor by x because in most instances r x for points near the axis. 

Similarly, from the solution of equation ( 30) we determine the formula 
for a constant continuous line source in a wind: 


Y = 

^ 4ffD 


/ 


[x-u(t-t* ) ] Vv^ 
‘ 4D(t-t») 


t-t» 


dt* 


Again, assuming the source is maintained indefinitely and making a change of 
variable to 


_ 1 

4D(t -t») * 


we have 



where Kg is Uie modified Bessel function of the second kind. 

For sufficiently large values of the argument pu/ 2D, we may use the 
asymptotic expansion for Kg, i.e. , 
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so that equation (44) becomes approximately 
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If <me is interested in the concentration relatively far downwind and near the x 
axis, which is usually the case, then a further simplification is possible. 
Expanding 


p = (x^ + /) 
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the as 3 m:iptotic solution becomes 
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“ 4Dx 
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In this expression one may also replace p by x without serious error. 

Thus far in treating continuous sources with or without wind, we have 
implied an atmosphere continuous in all directions, i. e. , one with no boundaries. 
Often, however, the source of the diffusing cloud is on or near the Earth’ s 
surface, where reflection, absorption, or deposition of particulate matter under 
gravity can occur. Although in reality all three effects take place more or less. 
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in many problems it can be assumed without serious error that the surface acts 
as an impervious boundary, implying total reflection. Since the continuous 
point source in a wind at or near the ground is of considerable importance, we 
will obtain the solution to this problem in the presence of an impervious boimd- 
ary ( the ground) . 

The method of images will be employed, as was done for the extended 
instantaneous plane source, hi obtaining the solution, equation ( 43) , the origin 
of coordinates was placed at the source. For reasons of convenience and 
S3Tnmetry we now place the origin on the ground and the source on the Z axis 
at a height z = H; therefore, the image source is placed at z = -H ( Fig. 6) . For 
the two sources, the solution is still given by equation ( 43) except now the 
distance to the observation point r from the actual source is given by 


r* = + y^ + ( z - H) ^ 



i 




and that of the image source is given by 




R* = x* + y^+(z + H)^ 


Taking the sum of the two solutions, we have 
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which is the exact solution. Making the same approximation for points near the 
X axis as was done previously, we have 

_ u(z+H)^ 

4Dx 4Dx 

e + e 
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For a source on the ground, one puts H = 0 in these expressions. Along the 
centerline of the plume, y = z = 0, this equation reduces to 


Q 4Dx 

D. Source Solutions in an Anisotropic /Vledium 

For a nonisotropic and m<ndng medium, the basic equation is equation 
(17) . Assuming that the velocity of the medium is parallel to the X direction, 
the equation in expanded form is 
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To obtain solutions of this equation appropriate to an infinite medium such as ! 

the atmosphere* ure flrst make a change of coordinates defined by ! 
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Then, by differentiation, the derivatives with respect to x transform as follows: 
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9z* 3z»^ D 3 


Also, the velocity along the new coordinate Is found by differentiating x’ with 
respect to time, resulting In 


n .^2 


Thus, the diffusion equation takes the form 


I2L + ui 
8 t 



(47) 


which has the same mathematical structure as the equation for isotropic media 
with an effective diffusivity equal to (D 1 D 2 D 3 ) Therefore, the appropriate 
solutions can be written by inspection in analogy to the solutions developed for 
isotropic media. 

The procedure is as follows. Treating the instantaneous point source 
first, we see that the solution of equation (47) is given by equation (42) ; 


X = 


Xo 


(27r)^'^V*’ 



(48) 



where now the standard deviation tr* , by analogy to equation ( 22 ), is defined 
by 
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We can write this in terms of the separate standard deviations for each coor- 
dinate direction, defined by the following quantities 


(T ^ = 2Dit 
X 


(T * = 2D,t 
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a ^ - 2D3t 
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Eliminating Dj, Dj, D 3 , we find 


o’ = (a (T (7 ) 

X y z 


V3 


(49) 


Also, eliminating Dj, Dj, D 3 from x’ , y’ , z’ , and u' , we have 


iCf = 


(a a cr ) 
X y z 
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(a cr (T ) 

. X y . 
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u» = 


(51) 
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(or 0 - cr ) 

X V z 


‘/3 


rinally, substituting equations (51), (50), and (49) into equation (48), we find 
the instantaneous point source solution for a moving anisotropic medium in terms 
of the original coordinates x, y, and z: 


Xo 

cr cr 
X y z 



Relative to an observer moving with the medium, u = 0; therefore. 


X = 


Xo 


(2ir) a o 
X y z 



(52) 


Similarly, for the instantaneous line source 



2 

27r(cr cr cr 
X y z' 



Since these solutions exhibit Gaussian distributions in the flow directions, they 

are basic forms employed in turbulent diffusion, with additional assumptions on 

O’ , (7 , and cr whenever Gaussian dispersion can be assumed. 

X y z 
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Many attempts have been made to solve the diffusion equation with 
variable diffusivities and D 3 . These efforts have succeeded only for 

certain special functional forms for these coefffeients. Details and references 
to the original investigations can be found in Reference 4. 
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IV. TURBULENT DIFFUSION 
A. Taylor's Theorem 

The basic solutions of the diffusion equation developed in the preceding 
sections are generally valid for molec\ilar diffusion. In general, however, 
dispersion of matter can also occur by turbulent motion. For example, the 
dispersion of a pollutant in the atmosphere almost always occurs in a turbulent 
flow field in which molecular diffusion is negligible. Unlike molecular motion, 
turbulent flow occurs on a macroscopic scale. Turbulence in the atmosphere 
can arise from thermal and pressure gradients, variations in wind velocity, 
boundary reflections such as occur on the grormd, and buoyancy forces which 
are related to vertical temperature variations. The complexity of turbulence 
is well known. Even with some simplifying assumptions, a theory of turbulent 
flow based on dynamical considerations leads to coupled differential equations 
which are virtually intractable to solution. 

Inasmuch as turbulence plays a primary role in airborne dispersion, it 

can be said flatly that the Fickian diffusion equations (17) and (ll), characterized 

by constant eddy diffusivities, fail completely in the atmosphere. The solutions 

of equation (17), the equation for a moving and anisotropic medium, have been 

employed as a basis to describe, with some success, the diffusion process in 

the real atmosphere. Even then, some modifications and some new assumptions 

are necessary for these solutions to conform adequately to empirical data. As 

we shall see later, a significant method of modification involves assumptions on 

the functional forms of c , o , and <r , and on the wind velocity, bi summary, 

X y z 

much of the empirical data can be absorbed in a theory that treats these param- 
eters not as fixed constants but as subject to variation in time or constant over 
limited regions of space. A fair question is how the solutions obtained from the 
molecular diffusion equation on the assumption of constant diffusivities can be 
carried over to turbulent diffusion in which the difffisivities are not constant. 
Obviously, if these are point functions of the coordinates, this cannot be done. 
However, if they depend solely on time, then as we have seen in Section n, a 
time scale transformation leads to an equation of the same mathematical form 
and, therefore, with the same type of solutions. 
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A fairly successful theory of turbulence has been based on stochastic 
theory rather than on dynamical considerations. A statistical technique which 
originated with Taylor and which has found frequent application is based on a 
Lagrangian treatment of the flow field, i. e. , one in which attention is focused 
on a single particle as it moves about the field. The theory involves the use of 
mean values of quantities characterizing the particle, such as its velocity and 
displacement, with fluctuations about the mean which are assumed to be sta- 
tionary, homogeneous, isotropic, and Gaussian. By stationary it is meant t iit 
statistical properties do not vary with time. For example, a covariance com- 
puted from the data in a particular time interval is assumed to be the same for 
any other time interval. Homogeneity implies that these properties do not 
depend on the space coordinates: that is, a measurement taken at a particular 
pcint in the field will r jssess statistical characteristics identical to those 
taken at some neighboring point. Isotropy means that the fluctuations do not 
vary in direction about a fixed point. 

With respect to the atmosphere, homogeneity in the horizontal direction 
is a reasonable assumption in regions under which the topography Is similar 
over a large area. Such cannot be said for the vertical direction. Because of 
the presence of vertical forces of gravity, buoyancy, and the Earth* s surface, 
the assumption of vertical homogeneity is unrealistic in most cases. The con- 
cept of isotropy suffers from similar limitations. 

Once a random process is assur.ed, a distrtbution function has to be 
postulated. It has been determined from empirical data that many atmospheric 
conditions can be represented by a Gaussian distribution, which yields the nor- 
mal probability curve. However, any other distribution may be assumed if it 
offers a better prediction of the variables under certain conditions. 


It is not the purpose of this report to give detailed accounts of the sta- 
tistical theories of turbulence. We limit ourselves to an elementary treatment 
of a significant statistical method founded on Taylor* s theorem and from it 
derive expressions for cr^, cr^, and cr^. This will be followed by two examples 

illustrating the use of these expressions in actual diffusion problems, one 
example from Brownian motion and the other from Sutton* s early work. Although 
this particular aspect of Sutton* s work has been superseded on theoretical 
grounds, it has produced several practical formulas with good predictive power 
under certain conditions. This section is ctmcluded with a brief description of 
the Hay-Pasquill method of cloud spread prediction. The purpose here is to 
provide a background for further reading. 
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We imagine an ensemble of particles in random motion and fix our atten- 
tion on a single particle. If the particle at time t - 0 is at the vector position 
T* and at time t at position r, the displacement T -T* will be a random function 
of time. If we now introduce a probability density function P (T -7* • t) , i. e. 
a probability per unit volume* then the probability that the particle is in the 
volume element dr about 7 at time t is 


P(7-7», t) d7 


The assumption of homogeneity implies that P is solely a function of the dis- 
placement 7 -7’ and not of the starting point* or release point* 7* . That is* 
every elementary region of the ensemble is represented by the same probability 
density. Physically* the meaning of P is that it approximates the fraction of 
particles in a given volume element. The term "approximates" is used because 
of the inherent statistical nature of P in particular and of the theory in general. 
If Q denotes the total mass released (which constitutes the ensemble) and x 
some mean concentration* it follows that 


x(7»* t) = QP(7-7»* t) 


( 54 ) 


If we place the observer in a frame of reference moving with the mean 
velocity of the ensemble* the mean velocity need not be considered. However* 
the mean of the square of the ensemble velocity is not necessarily zero* although 
it must be independent of time. Thus* if u denotes the x component of velocity 
and the bar over a quantity denotes its mean value* then 


u(t) = 0 u^(t) = constant 


Wc have similar expressions for the y and z components. After a small time 
T* the velocity will be u(t + t). Therefore* the velocity covariance is* by 
definition* 





u(t) u(t + t) 
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Since the process is stationary* the covariance is a constant. This implies that 
the starting point t at which u(t) and u(t r) are computed does not affect the 
mean value of their product. An important result of this is that the covariance 
is an even function of r* meaning that it has the same value for t as for -r. To 
understand this* let us take the starting time t « 0* so that the covariance is 


u(o) u(t) 


Now let us take t = -r* giving 


u(-t) u(0) 


Because of the assumption of stationarity these must be equal; 


u(0) u(t) = u(-t) u(0) 


which is clearly even. 

The next step is to introduce a velocity correlation function R( r) ( some- 
times called a Lagrangian correlation coefficient) deflned by the ratio of the 
covariance to the mean of the squared velocity: 


R(t) = 


u(t) u(t + t) 


u* 


( 55 ) 


which is also independent of t because numerator and demonator are Independent 
of t. 


We now compute the mean time rate of change of the square of the dis- 
placement x(t) : 
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Since 


u(t) , x(t) = f u(t') dt’ ; 
0 


hence, the above becomes 

t 

— x^(t) =2 f u(t) u(t‘) dt* 
6 


where we used the justifiable process of interchanging the order of integration 
and averaging. Making a change of variable defined by t* = t + t, we have 


^ x^(t) = 2 f u(t) u(t + t) dr 


t 

= 2 f u(t) u(t + t ) dr 
0 


where in this last step we replaced the integration from -t to 0 by 0 to t because, 
as noted previously, the integrand is an even function. Finally, substituting 
from equation ( 55) , we obtain 


I X*(.) = 


/ 


0 


R( t) dr 


which upon integration over the interval t results in 




_ t t» 

x*(t) = 2u* f f R(t) drdt’ 
0 0 


_ t 


= 2v? f (t - t) R(t) dr 


( 56 ) 


This result is known as Taylor’ s theorem and is of fundamental importance 
in turbulence. It relates the mean squared displacement of a particle in time t 
to its mean squared eddy velocity and velocity correlation function. 

We now seek a connection between the autocorrelation function R( t) and 

the standard deviations o , a , and <t which quantities, as we have noted, 

X y z 

characterize the spread of the dispersion along the three coordinate directions. 
We recall that in one dimension a was defined by 
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CO 

—CO 


Extending this definition to three dimensions, we have 


00 GO CD 

^ I f f (x, y, Z, t) dx dy dz 

X Q ^ 

-00 -00 -CO 


with similar expressions for o ^ and cr Substituting from equation (54), we 
get ^ ^ 


CO 00 CO 

ff f f f x^P(x, y, z, t) dx dy dz 

_oo -00 -CO 
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but from statistical theory this integral is just the mean or expected value of 
x^. Thus 



_ t 

= 2 u* f (t - t) R(t) dr (57) 

0 


from Taylor* s theorem. This is the relation sought. L(i general, for all three 
dimensions, we have 

_ t 

I ■ I") 

y 0 ^ 

_ t 

(T ‘ = 2 f (t - t) R (t) dr 
z ^ z 


where v and w are the velocity components in the y and z directions, and R^ 

and R are their respective correlation functions. Thus, in this theory the 
z 

problem of turbulence is reduced to that of determining the velocity correlation 
function. The difficulty, however, is that there is no satisfactory theory from 
which this function can be predicted. The function can be predicted if it can be 
assumed that the velocity history is a Markov process. A Markov process is 
one in which the value of a random variable at times greater than some t 
depends on the value it has at t but not on any previous value. There are some 
theoretical inconsistencies in the assumption of a Markov process for the 
velocity fluctuations, however. These will be discussed in a later paragraph. 

Some general properties of R( t) are worth noting. We have already 
seen that it is independent of t and an even function of t. Also, at t = 0, 

R(0) = 1, and as T — 0, R(t) — 0. In general it can be shown that for all t, 

-1 ^ R(t) ^ 1 
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B. Example from Brownian Motion 
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It Is an enlightening exercise to apply Taylor* s theorem to Birownlan 
motion, which Is the agent by wJiich molecular diffusion occurs, and to show that 
it gives the same result for a as does the phenomenologically based diffusion 
equation. The exercise Is expressive of the significance and power of the sta- 
tistical approach. When small particles of colloidal size suspended In a fluid 
are observed under a microscope, the particles are seen to move about with 
apparently random motion. This motion was observed by the botanist Robert 
Brown In 1826. Not until 1905, when Einstein published his classic paper, was 
the phenomenon completely understood. He showed that the motion la the result 
of collisions with the molecules of the surrounding fluid, which eventually leads 
to the dispersal of the particles. 

We *q)ply Newton’ s second law of motion to a single particle of mass m 
and velocity u moving In one dimension. The assumed forces acting on the 
particle are viscous drag and collision forces taken as a group of random 
impulses. The drag force Is assumed proportional to the first power of the 
velocity, with K denoting the proportionality constant and f(t) denoting the 
random force component along X. The equation of motion Is 


m ~ = -Ku + f(t) 
dt 


Writing ^ = K/m, we have 



m 


This first-order, linear differential equation can be integrated with the aid of 
the integration factor e^. Multiplying by this factor and rearranging gives 


dt 


(ue 


flii g/3t 

' m 


* 


55 



r 


I 


i 


I 



I 


k 


so that 






dt» 


where Uq is the initial velocity. Solving for u. 


-/3t 

u = U(,e 


+ 



/ f(t») e^’ dt» 
0 


(58) 


This result shows that the particle* s velocity depends on the initial velocity and 
on the random acceleration f/m. For times t » 1/ 0, the first term decays to 
zero, making u independent of its initial velocity. It may be stated that after a 
sufficiently long time the particle ’’forgets** its initial velocity and is thereafter 
’’pushed” around by the random forces. Thus, after a long time the process 
becomes a Markov process. It might also be noted that in spite of the presence 

of the decay factor e ^ in the second term, this term does not necessarily 
decay, at least at the same rate, as the first term because the integral is not 
constant but some function of t. 

If we now take the ensemble average of equation (58) (the average of a 
large number of randomly moving particles) , the contribution to u from the 
second term vanishes; therefore. 


-0t 

U = Uq e 


(59) 


Moreover, from the principle of equi -partition of energy founded on the statis- 
tical theory of gases, we have 
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where the left side is the mean kinetic energy of the particle, k is the Boltzmann 
constant, and T is the absolute temperature. This simple result shows that 

u^ is independent of time, indicating that Brownian motion involves a stationary 
process. 

{ 

We are now able to determine the velocity correlation function for 
Brownian motion. Because the process is stationary, we can take t = 0. Thus, 

I . with the aid of equation (55) , we find 




Substituting this correlation function into Taylor* s theorem, we have 
— t 

(T^ = 2 u* f (t - t) R(t) dr 






> 




To interpret this resiilt and relate it to the expression obtained from the phe- 
nomenological theory, we need some insight as to the magnitude and significance 
of the quantity /3. Its reciprocal has the dimensions of time. It is the time for 
the velocity to decay to 1/e of Its initial value, or the "relaxation time" of the 
viscous effects. Its magnitude may be computed from Stokes' law: 

^ m 


where a is the radius of the particle and /i is the viscosity of the fluid. For a 
particle of mass m =* 10“^ g In air, = 1.6 x io“* g-sec“*-cm“\ and a lO"* 
cm, we have 

/3“*« 1.4 X 10"* sec 

which is much smaller than the time scale of normal diffusion. Thus, for 
t » 1//3 we neglect the second term in equation (60), giving 


a 


2 


X 



t 


Comparing with cr * = 2Dt, we see that both expressions are in accord in that each 

X 

depends linearly on t. Also, it follows that 



JcT 

Sirtifi 


This relation has been repeatedly verified by experiment and is a well established 
result of the theory. 
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We have seen that for Brownian motion the velocity correlation function 
has the form of an exponential decay and, therefore, the motion is of the Markov 
type. If we assume that eddy motion is also of this type, then 


R(t) ~ e 


-t/ 0 


(61) 


where 0 would be some appropriate time scale for the process. For Brownian 
motion, 0 is detemined by viscous effects which was found to be 0 = 1/^. It 
would seem, therefore, that 6 in general is determined by some characteristic 
length related to the size of an eddy and some characteristic velocity. On this 
basis, 0 would be expected to be proportional to the ratio of the characteristic 
length to the characteristic velocity. All this appears quite reasonable; how- 
ever, an obvious difficulty arises from the assumed functional form of R(t) if 
we recall that it should be even. Replacing t with its absolute value in equation 
(61) makes it an even fimction; however, it has a discontinuity at t = 0. Physi- 
cally, this implies an infinite acceleration of the particle (or eddy) or a vanish- 
ing inertia. Thus, the type of function being considered is not acceptable on 
theoretical grounds. However, for large values of t the particle is no longer 
influenced by its initial velocity; therefore, a Markov process should be an 
acceptable model for such times. In the limit of large t, the dispersion <r^ is 
proportional to the time, which is characteristic of molecular diffusion but not 
generally valid in turbulence. The obvious inference is that the Markov process 
is applicable only when the turbtUent medium has reached a quiescent state 
dominated by molecular diffusion, at which time the problem may no longer be 
of interest. 

As indicated by Sutton 15) , the difficulty with an exponential decay 
correlation function is that it implies eddies all having the same size, analogous 
to the molecules of a simple gas, which is definitely not in accord with the facta 
of turbulence. The conclusion is that this type of function cannot form the basis 
of any general theory. 

Other theories, including the "mixing length," have been proposed. In 
these theories ouw introduces a characteristic length, akin to the mean free path 
of a molecule whose magnitude depends on the intensity of turbulence. The idea 
is that a randomly moving fluid element transfers its momentum to the mean 
flow in a distance of the order of the mixing length. This implies that the ele- 
ments behave independently of the mean motion for brief periods during which 
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they transfer their momentum, much like independent, colliding molecules. No 
model of the structure of an eddy is necessary. The velocity fluctuation there- 
fore becomes a function of the mixing length. Unlike the mean free path, the 
mixing length may be a function of position, mean velocity, and other variables. 
In practice, the mixing length is a convenient parameter lacking any solid physi- 
cal foundation. A result of the theory is an expression for a virtual coefficient 
of diffusion equal to the product of the mixing length and the square root of the 
mean of the squared velocity fluctuation. The details of the derivation can be 
found on page 72 of Reference 5. 


C. Sutton's Formulas for the Dispersion 

In Sutton* s analysis, it is assumed that the correlation between fluctu- 
ating velocities near a sir ooth surface depends on the mean eddy energy p u*, 
the fluid viscosity p, and the time t. The only dimensionless number that can be 
formed with these quantities is p/p i? t, which, with v = p/p, can be written 
v/ u2 T. The quantity v is the kinematical viscosity. Considering the limiting 
properties of R(t) for small and large r, the simplest function with these prop- 
erties is the power law 



where n is a positive number. Thus, in nonisotropic diffusion, e. g. , diffusion 
near the ground, one can assume 
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where u, v, and w are the velocity fluctuation components. Furthermore, 
Sutton suggests that for rough surfaces v In these expressions be replaced by 
the ’ 'macroviscosity,” a quantity without precise definition except that It is, 
for the most part, empirically determined and several orders of magnitude 
greater than v. Its magnitude, from Sutton' s computations, ranges from 0.016 
for ice and mud flats to 560 for thick grass up to 50 cm high. 

From equation (57), we calculate 

a ^ = 2 u^ f ( t - t) / = — \ dr 

^ 0 V + u^t/ 

I t t 

.nr dr nr t dr 

tv I n - V I II ' 

6 ( + u^ t) ” 6 ( + u^ t) ” ^ 

_ v + u^ t)^ ^ 2i^ 2vt 

(l-n)(2-n)u^ (l-n)(2-n)u^ ^ ” 


A 

Neglecting terms of order v compared to u t, we have 


a 2 = 


X 


„ n - 1 ^\2-n 

t) ^ 

(1 - n)(2 - n) u^ 


Similarly, we find 



2v (v^ t) 

(l - n)(2 - n) v^ 


cr = 
z 


.. n ,\2-n 
2 _ 2v (w^ t) 


(1 - n)(2 - n) w^ 
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Defining ’’generalized diffusion coefficients” C^, C^, by 


C = 


C = 

y 



c 

z 


4p 


(1 -n)(2 -n) w2 



where U Is the mean (constant) cloud velocity, one can write 


2 

(T 

X 








1 2 / \ 2^n 

(Ut) 


(62) 


For n » 1, one finds that Is proportional to t, a behavior characteristic 
of molecular or Brownian diffusion. For n = 1/4, a* Is proportional to 
which, according to Sutton, gives a reasonably accurate description of diffusion 
In the atmosphere from a few meters to hundreds of kilometers. 

We can arrive at Sutton’ s formula for the Instantaneous point source In 
a reference frame moving with the cloud by formally substituting the relations 
of equation (62) Into equation (52) : 
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From this basic solution Sutton obtains, by integration, the formula for the 
continuous point source near the ground (the ground reflection is taken into 
account) : 



y z 


where x = Ut. Expressions for line sources are also given (see page 288 of 
Reference 5). 

A first criticism of Sutton* s development is that the constants C , C , 

X y 

C do uot have fixed dimensions but depend on the value of n. Furthermore, the 
z 

correlation function varies as t , which yields divergent diffusion time scales 
and Infinite spectral density at zero frequency under Fourier analysis. In spite 
of these theoretical inconsistencies, Sutton’ s basic formula, equation ( 63) , is 
found to be identical to that derived by some other theoretical methods, e. g. , 
by treating turbulent diffusion as a random walk process. Csanady [GJ suggests 
that the theoretical discrepancies may be due to the nonuniformity of the tur- 
bulent field near the ground as the cloud rises and the eddies grow in size. 

Since the field is no longer stationary, the concept of a velocity correlation 
function is undermined. We might reverse the reasoning as follows. Starting 
with the presumably correct power law formulas, equation (62), a velocity 
correlation function deduced from tnem would be invalid. 


For a continuous elevated point source in the presence of a constant mean 
wind along x and total reflection at the ground, Sutton assumes the expression 




^ I 

1 I ) 




Y = ■ * — — 

^ 27rUa (T 

y 2 


where H is the source height measured from an origin on the ground beneath the 

source. !h this formula we have replaced Sutton* s C and C with a and a 

y z y z 

defined by equation (62) . This general form for the concentration allows differ- 
ent assumptions on the values of and and has been widely employed by 

others to predict cloud concentration at ground level from smoke stacks. At 
ground level ( z = 0) , the formula becomes 




2cr * 2 <t 

\ y 2 


and along the axis of the plume (y ~ O) , it becomes 


Y = e 

* ttUo a 


It is easily seen that the maximum concentration along the plume axis occurs 
at the point downwind where the relation a - H/^/T satisfied; the concentration 
there is 


2Q / i 
^max irUerf I a 
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D. Hay*Pasquiil Method for Determining the Dispersion 

In deriving Taylor* s theorem it was explicitly assumed that the velocities 
were Lagrangian velocities, i. e. , the mean velocity and fluctuating components 
of a particle as it moves about from one point to another. In contrast to this 
Lagrangian description, the Euler description deals with velocities at a fixed 
point of the fluid. Since in the determination of experimental data on diffusion 
velocities samples are often taken with instruments fixed at a point in the flow, 
it wo old be useful to determine the relationship, if any, between Lagrangian and 
Euler velocities. The difference is, primarily, that the velocity fluctuations at 
a fixed point occur in a shorter time than those of a drifting particle. The Hay- 
Pasquill method is based on the realization of this time difference. 

In this method one introduces a shortened time interval T, or simulation 
time, related to t by t = /?T, where p is some empirically determined constant. 
With this new time scale and by the same arguments used to derive Taylor’ s 
theorem, one now obtains 





T t’ 

f f R_( t) dr dt’ 
6 6 


for the simulated mean-square displacement, where U is the Eulorian velocity 

£ 

and R ( t) is the Eulerian velocity correlation function. One then assumes that 
£ 

the Langrangian velocity fluctuations arc equal to the Eulerian fl uctua tions but 

that they occur on a time scale Ic^fror by a factor of p. That is, U_^ = U*, 

_ £ 

where U* is the corresponding Lagrangian quantity. Employing this equa^y in 

equation ( 56) and then equating equation ( 5f>) to the above expression for x^ , 

it follows that 


Rg(r) = R(Pt) 


where the right member is the Lagrangian correlation function. This result 
implies that the Eulerian and Lagrangian correlation functions have the same 


{ 


( 


I 



t 


65 


1 


1 


I 


shape but that the Lagrangian one falls mote slowly. Although it. has been 
definitely determined that the two correlation functions do not. In reality, have 
the same shape, the standard deviations obtained by this method give good 
predictions of cloud concentration [ 4, 6] . 
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V. NASA/MSFC MULTILAYER DIFFUSION MODEL 




A. Point Source Solutions in a Layer 
with Impervious Boundaries 


The real atmosphere quite often possesses a stratified structure, pri- 
marily as a result of temperature and wind shear. The turbulent mixing of a 
cloud with ambient air (entrainment) tends to decrease the buoyancy and vertical 
momentum of the cloud. As a result the region in which turbulent mixing is 
occurring will be bounded by a stable layer which acts as a lid, preventing com- 
plete admixture of the air across the boundary. This would be the case when 
buoyancy foices are stabilized by gravity. Given that the ground also acts as a 
similar boundary, the task is to determine solutions in a vertical layer of 
atmosphere. The simplest diffusion model w'hich accounts for such boundaries 
is one which assumes total reflection at the boundaries. The solutions for the 
instantaneous and continuous point sources will now be derivev-*, followed by a 
summary of the multilayer model developed by the Cramer Company for NASA. 

A layer, by definition, consists of two distinct boundaries. Mathemati- 
cally, the problem can be treated by the method of images. The flux reflected 
at the first boundary is reflected at the second boundary which, in turn, is 
reflected at the first boundary, etc. The general solution can be constructed by 
assuming multiple images, one for each reflection. Each image upon reflection 
at the opposite boundary creates a second image which, in turn, creates a third, 
etc. The essential steps will be given in the following paragraphs. 

The origin of coordinates is placed on the ground, z = 0, the source at 
z = H above ground, and the top of the mixing layer at z = h, with h > H. Total 
reflection is assumed at both bouridaries, the ground and the top of the mixing 
layer. For the Instantaneous and continuous point sources we assume the basic 
Gaussian form 




2(7 


X = Ae 
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where we let 


A = 


Q 

a O’ 

X y z 



for the instantaneous source* and 


A = 


Q 

2irUo O’ 

y 2 



e 


y 


for the continuous source. 

Consider first the reflection of the real source at the boundary z= 0. For 
this reflection an image source is placed at z = >H. The solution due to this 
image is 



The flux from this image is also reflected at the boundary z - h. Therefore, 
we place a second image at z = 2h + H (at the same distance from the boundary 
z = h as the first image) . The solution due to this second image is 

fz-(2bfH)l^ 

2<t^ 



The flux from this secondary image is itself reflected at the boundary z = 0; 
therefore, a third image is placed at z = -( 2h + H) , giving rise to the solution 



The flux from this tertiary image is reflected at z = h; therefore, a fourth image 
is placed at z = 4h -i- H, giving the solution 



A fifth image at z = -( 4h + H) yields 



Generalizing the results and taking the sum of the solutions, we have the expression 




To this expression we now have to add the sum of the solutions arising from 
reflection of the real source with the boundary z » h. In this case we first 
place an image source at z » 2h - H, yielding the solution 



For reflection of this image at z = 0, we place a second image at z = -(2h - H) , 
yielding the solution 



/ 

For the next two reflections at z = h and z = 0, we have, respectively, 




etc. The sum of these solutions is 



A Z « ^ * Z • 

1=1 1=1 
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The term i = 0 is not included here because it has been included in the preceding 
partial sum, equation (64). 

Thus, the general solution is the sum of the two partial sums, equations 
(64) and (65): 



In the literature on diffusion these sums are written in several equivalent forms. 
The form used by the Cramer Company in its reports (7, 8J is obtained by 
writing the terms for I = 0 explicitly and combining the two sums into one as 
follows: 




Except for minor dlfferencest such as the use of the equivalent (2ih -i- H - z)^ 

for (z - H - 2ih)^, and the use of H for our h, the expression in braces is 

m 

identical to the ’’vertical term” of the Cramer group. 

Another equivalent form of equation (66), and the most concise, is 
obtained by replacing the summation range in the second sum by i » -1 to 
and accommodating this change by an appropriate sign change in the affected 
terms. One can easily verify the result 


CO 


X = A Yj 

i=-OD 


(z+H->-2ih)^ (z-H+2ih)^ 

2(t * 2(7 ^ 

z z 


( 68 ) 


The ground-level concentration along the axis of the plume is of interest. 
With z = 0 and y = 0, we obtain in the case of the continuous source: 




_ (H-f2ih)' 


CD 


2(7 


^ " ffU(7 (7 ^ ® 

y z 


(69) 


We recall that characterizes the spread of the concentration in the 

vertical direction and, therefore, increases with time; for example, as t -» <>o, 

(7 00 , and X approaches some asymptotic vali» independent of z. Since the 

z 

total amount of diffusing matter must be conserved, the integration of equation 
( 68) over all space must equal the total mass released. The integration results 
in the asymptotic form 
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Thus, the summation factor Ir. equation ( 68) , which expresses the vertical 

distribution, approaches the distribution 2tt or / h asymptotically. That is, 

z 

the Gaussian distribution eventually becomes rectangular. 


B. Summary of the NASA/AASFC Model 

The NASA/ MSFC Multilayer Diffusion Model developed by the Cramer 
Company is employed in the predietion of fuel hazards resulting from NASA 
activities. In particular, this model was designed to permit concentration and 
dosage calculations downwind of toxic clouds from rocket vehicles. The basic 
concepts of this model are summarized here in the context of the stated purpose 
of this report. The specific details of the method and associated computer 
program can be found in References 4 and 6. 

The stratification of the atmosphere into regions of significantly different 
meteorological parameters, such as wind velocity and temperature, is the basis 
for assuming a multilayer model. Each layer is assumed homogeneous and the 
boundaries impervious to the turbulent flux. However, the formulation does 
permit taking into account flux of matter across the boundaries due to gravita- 
tional settling and precipitation scavenging. The diffusion process is repeated 
from layer to layer, each layer assigned a new set of source and meteorological 
data. 

The model takes into account the loss of matter from the cloud resulting 

from decay processes, precipitation scavenging, and gravitational settling. 

All these effects combine to deplete the matter that forms the cloud. The 

expression for the concentration, therefore, consists of the product of two 

terms: a diffusion term and a depletion term. If a Gaussian distribution is 

assiuned, the diffusion term is that developed here, equation (67), for the 

instantaneous or continuous point source, depending on the function A as already 

defined. The standard deviations a , a , and a are now ganerallzed entities 

X y z 

significantly, but not conceptually, different from the relatively simple meaning 
we have attached to them. They are defined by certain rather awkward seml- 
cmpirical relationships involving quantities that vary during the time required 
for the cloud to stabilize with respect to ambient air. In the main, therefore, 
these quantities are functions of one or more of the following: ( l) the standard 
deviations of azimuth and elevation angles of the wind; (2) the distances over 
whieh vertical and crosswind expansion of the cloud occurs; (3) the standard 
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deviations of the vertical, crosswind, and alongwind distributions at the time 
of cloud stabilization; ( 4) vertical and lateral diffusion coefficients; and ( 5) the 
time for the cloud to reach equilibrium with ambient conditions. Much of the 
input data required to calculate these cr^ s must b^ empirically determined. 

If gravitational settling is neglected, the depletion term consists of the 

-bt 

product of two exponential decay factors of the form e » one for decay processes 
and one for precipitation scavenging, where b is some appropriate constant for 
the particular process. 

If gravitational settling is not neglected, then the sum (but not A) in the 
diffusion term, equation (68), is replaced with the expression 

’ 2o"‘5 

Z . Z 


where V is the settling velocity. This replacement is made because the ground 
can no longer be considered a reflecting boundary, reflection now occurring 
only at the upper bovindary z = h. Equation (71) is, therefore, really a diffusion 
term. Neglecting for the moment the term Vx/U, equation ( 7l) can be written 
as 



2cr “ 
z 

+ e 


(72) 


which, aside from the factor A, we recognize as the solution for the point source 
with source at z = H and a single reflecting boundary at z h. 

The inclusion of the term Vx/ U is based on Schmidt* s studies on sedi- 
mentation as modified by Csanady 1 6] . In equation ( 72) the source height H is 
replaced by H - Vx/ U, resulting in equation ( 71) . 
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